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l' Abstract 

Q ' We solve a weakly singular integral equation by Laplace transforma- 

S—j \ tion over a finite interval of R. The equation is transformed into a Cauchy 

CZ2 . integral equation, whose resolution amounts to solving two Fredholm in- 

tegral equations of the second kind with regular kernels. This classical 
^ , scheme is used to clarify the emergence of the auxiliary functions express- 

ing the solution of the problem. There are four such functions, two of 
y\i ' them being classical ones. This problem is encountered while studying 

the propagation of light in strongly scattering media such as stellar atmo- 
spheres. 

Key words: Weakly singular integral equation, finite Laplace trans- 
form, sectionally analytic function, radiation transfer theory, stellar at- 
mospheres. 

1 Introduction 

The integral form of the equation describing the radiative transfer of en- 
ergy in a static, plane-parallel stellar atmosphere is pQ 



[ K(T-T')S(a,b,T')dT', 
Jo 



S(a,b,T) = S (a y b,r) + a K(r - T ')S(a,b,T')dT' , (1) 



1 



where S is the source function of the radiation field and So describes the 
radiation of the primary (internal or external) sources. These functions 
depend on the two parameters of the problem: the albedo a g]0, 1[, which 
characterizes the scattering properties of the stellar plasma, and the op- 
tical thickness b > of the atmosphere. They also depend on the optical 
depth t 6 [0, 6], which is the space variable. Equation Q means that the 
radiation field at level r is the sum of the direct field from the primary 
sources, and the diffuse field having scattered at least once. 
In the simplest scattering process conceivable - i.e., a monochromatic and 
isotropic one - the kernel of the integral equation (0 is the function 

K{j):=\M\r\) (reR*), (2) 
where E\ is the first exponential integral function 

Bi(r):= [ exp(-r/x) — (r > 0). (3) 
Jo x 

Since _Bi(jr|) ~ — ln(|r|) when |r| — *• + , the kernel of the integral equa- 
tion Q is weakly singular on its diagonal. The free term So includes the 
thermal emission of the stellar plasma - of the form (1 — a)B*(r), where 
B* is a known function - and the contribution aJo xt (b,r) of the external 
sources via the boundary conditions: see the introduction of PQ. Hence 
So{a,b,r) = (1 — cl)B*(t) + aJ^ xt {b,T). In a homogeneous and isother- 
mal atmosphere assumed to be in local thermodynamic equilibrium, the 
function B* coincides with the Planck function B(T) at the (constant) 
temperature T. Moreover, J§ — in the absence of external sources 
and thus 

S (a,b,r) = (l-a)B(T), (4) 

which shows that So is independent of r in this model. The solution S to 
the problem Q is then 

S(a,b,T) = (l-a)B(T)Q(a,b,r), (5) 
where Q solves the following integral equation: 

Q{a,b,T) = l + a I K(t - T')Q(a,b,T')dT '. (6) 
Jo 

It is proved in [5] that the space C ([0, &]) of the continuous functions 
from [0, b] to M is invariant under the operator 

A : / - A/(t):= F K(t - r')/(r')dr', (7) 
Jo 

with norm 

fb/2 

||A||oc=/ Ei(r)dT=l-E 2 (b/2). (8) 
Jo 

Here, ^(r) := L exp(— r/x)dx is the exponential integral function of 
order 2. Equation @, which can be written in the form 

Q = l + aAQ, (9) 
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has therefore a unique solution in C°([0, b]) provided that a < 1 or 
b < +00. 

This problem is a basic one in stellar atmospheres theory, and more gen- 
erally in transport theory. It describes in the simplest way the multiple 
scattering of some type of particles (here photons) on scattering centers 
distributed uniformly in a slab of finite thickness, a very simple ID- 
configuration. Its applications in astrophysics and neutronics - among 
other fields - are presented in 0. It is important to solve this problem 
very accurately, thus providing a benchmark to validate the numerical 
solutions of integral equations of the form Q. 

Physicists and astrophysicists have developed many methods for solving 
integral equations of the form (0 with a convolution kernel defined by (J5| 
The main steps for solving the prototype equation JSJ are summa- 
rized in a recent article |5j, which contains accurate tables of the function 
(1 — a)Q for different values of the parameters a and b. While reading this 
paper, one is struck by the complexity of the "classical" solution to Eq. 
iJSJl. It requires the introduction of many intricate auxiliary functions in- 
troduced in the literature over more than thirty years. The reader quickly 
loses the thread of the solution, which reduces his chances of exploiting it 
for solving problems of the more general form Q. 

The aim of the present article is to get around this difficulty by solving Eq. 
© straightforwardly, introducing as few auxiliary functions as possible to 
express its solution. These functions are briefly studied in Appendix A. 
The method is based on the finite Laplace transform, which reduces the 
problem © to solving a Cauchy integral equation over [— 1, +1], which in 
turn can be transformed into two Fredholm integral equations over [0, 1]. 
This approach has been developed in transport theory after the publica- 
tion in 1953 of the first English translation of Muskhelishvili's monograph 
Singular integral equations see, e.g., |7j, [Sj and It can be consid- 
ered as an extension of the Wiener-Hopf method |10l for solving integral 
equations of the form Q with b < 00. Both methods are characterized 
by an intensive use of the theory of (sectionally) analytic functions, which 
allows solution of Eq. Q in a concise manner. This is obvious when 
comparing the solution we derive here to the classical solution of the par- 
ticular problem J5J, which does not use any calculation in the complex 
plane. The former method clarifies the origin and the role of the auxiliary 
functions expressing the solution to a problem of the form Q . Since these 
functions are independent of the source term 5*o, they are "universal" for 
a given scattering kernel. 

The remainder of this article is organized as follows: in Sec. 2, the finite 
Laplace transform of the Q-function is calculated on the basis of some 
recent developments on Cauchy integral equations |11 | -|13 | . Then the 
Laplace transform is inverted and the solution achieved in is concisely 
retrieved with the help of the theorem of residues (Sec. 3). It involves 
two functions F+ and F- with remarkable properties, as shown in Ap- 
pendix B. The difficulties arising from the numerical evaluation of the 
latter functions are investigated in |S]. 
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2 The calculation of the finite Laplace 
transform of the Q-function 

Supposing < b < oo, we plan to solve Eq. © by Laplace transform 
(LT) over [0,6]. This operator is denned on C° ([□,&]) by 

£/(«):= / /(r)exp(-Tz)dr (z G C). (10) 
Jo 

Since b < 00, the finite LT of a continuous function is defined and analytic 
in the whole complex plane. The inversion formula 

f(r)=^{ £/(z) exp(rz)dz (11) 

is valid at any r G]0, &[, with no restriction on c £ R. The symbol ^ 
on the right-hand side of Eq. mill means that the integral is a Cauchy 

c + iX 

r c+ioc 

principal value at infinity, i.e., + — lim 

c-iX 

With the intention of taking the LT of both members of Eq. we 
note that the LT of the function K, as defined by Eqs. exists on 

C \ {±1} and is 



^)=-(lA)-i/ 1 I ^-exp(-6 Z )^ 1 
where w : C \ {±1} — * C denotes the function 



expt-6/.)^, (12) 



The three integrals on the right-hand side of Eqs. 1121 - 1131 are Cauchy 
principal values over ]1, +00 [, ] — 00, — 1[ and ] — 1, +1[, respectively. From 
the definition Q of the A-operator, we infer that the finite LT of A/ is 

£(A/)(*) = w(l/z)£f(z) -if £/(!/«) "'" 



2 Jo 1 — zu 

-exp(-fe)i f £/(-l/«) exp(-6/u)-£U (14) 
2 Jo 1 + z« 

In addition, the finite LT of the unit function is z — > (l/,z)[l — exp(— bz)]. 
Taking the LT of both members of Eq. © and changing z into 1/z, we 
obtain the following integral equation for £Q: 

rl j 

T(a,z)£Q(a,6, 1/2) - / £Q(a, 6, 1/w) — — 
2 Jo U — 2 

a f 1 du 

+ exp(—b/z)—z / £Q(a, b, — 1/u) exp(— blu) = co(z), (15) 

2 Jo u + z 

where 

T(a,z):=l-aw(z) = l-±z (16) 
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and 

co(z):=z[l-exp(-b/z)]. (17) 

The function T is our first basic auxiliary function. Its main properties are 
summarized in Appendix A. The function Co satisfies the three following 
conditions: (t) it is defined and analytic in C*, (ii) it is bounded at infinity 
(with limit equal to 6), and (Hi) it satisfies, in a neighborhood of 0, the 
conditions 

lim co(z) < oo , lim [cq(— z) exp(— b/z)] < oo. (18) 

z—*0 z — >0 

Re(z)>0 Re(z)>0 

Replacing u by — u in the second integral of Eq. H15I . we obtain a Cauchy 
integral equation on [— 1, +1] for the function z — > £Q(a, b, 1/z). Since the 
latter is not Holder-continuous at 0, it seems inappropriate to undertake 
the resolution of this equation. A better approach consists in observing 
that the functions z — > £Q(a, b,l/z) and z — > £Q(a, b, —1/z) x exp(— b/z) 
are solutions of two coupled Cauchy integral equations on [0,1], which 
can be uncoupled by adding and subtracting them. The resulting Cauchy 
integral equations are then reduced to two Fredholm integral equations 
of the second kind on [0, 1] with regular kernels. This general approach 
was first introduced in transfer theory by Busbridge |7] and developed by 
Mullikin et al. 8 9 and Rutily et al. 13 . The last mentioned reference 
gives a synthesis that the reader may consult for details. It contains the 
proof of the following result, which we omit here: the unique solution, 
analytic in C*, to an integral equation of the form 1151 with free term 
satisfying the conditions (i)-(iii), can be written in the form 

£Q(a,b, l/z) = -[u-(a,b, z)rjo,+ (a,b, z) + u+(a, b, z)rjo,-(a, b, z)], (19) 

where, for any z £ C \ iR, 

it± (a, b,z):= Y[^(z)]H(a, z)£±(a, b, -z) 

T Y[-$t(z)]H(a, -z)C±(a, b, z) exp(-b/z), (20) 

and 

1 z /' +io0 , 
Tj ,±(a,b,z) := -co, T (a,6,+0) + — / H(a,l/z) 

2 2l7T J_ ioa 

x (( ± )-(a,b,-l/z')c (-l/z')-^. (21) 

1 + zz 

Equation 1201 involves the Heaviside function Y, which vanishes over R* 
and is equal to unity over R+. On the right-hand side of Eq. 121H . we 
have introduced 

c , T (a,b,+0) := lim [c (z)±c {-z)cxp(-b/z)]. (22) 

Re(z)>0 
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Since this limit vanishes for the particular co(z) given by 1171 . Eq. 12111 
can be rewritten as 



r) ,±(a,b,z) = -— J H{a, l/z')((±)~ (a, 6, -l/z) 

^-^')]^ W y (23) 

The functions H(a,z), <^+(a,b,z) and C,-{a,b,z) appearing in Eqs. 1201 
and 1231 are the other three auxiliary functions of our approach. They 
are defined in Appendix A, which synthesizes their main properties. The 
//-function is a classical function of radiation transfer theory, recently 
reviewed in |12|. Most results concerning the functions (± are summarized 
in |13|. They are sectionally analytic in the complex plane cut along the 
imaginary axis, with limits (C±) + an d (C±)~ on the right- and left-hand 
sides of this axis respectively. These limits satisfy the relation 

H(a,z )((±)~(a,b, -zo) = ^H(a,-z.o)((±y(a,b,zo)exp(-b/zo), (24) 

for any zo £ iR*. This can be seen by putting z — > zq on the right, 
then on the left, into Eq. 1201 . taking into account the continuity of the 
functions u± on iM |13|. The limits of the functions £± on both sides of 
the imaginary axis also satisfy the relation 

H(a,z )((;±) + (a,b, -z ) = H(a, Zo)((±)~(a, b, -zo) 

± H(a, -2o)(C±)+(a, b, zo) exp(-b/z ), (25) 

which follows from 1A111 and Plemelj's formulae 6 . 

In the expression 1231 of the functions 770, ±, replace the term containing 
exp(bz') by its expression given by 1241 (with zq — —l/z'). One obtains 



-too 



??o,±(a,M) = — 77T- I H{a,l/z')((±) (a, 6,-1/2') 

2 "T J-ioo 

^VTTTz ± V 1 ni^- ^ 

This integral may be calculated using the contour of Fig. 1. The integral 
along the half-circle T e tends to — (l/2)H(a, oo)£±(a, 6, — oo)(z =p z) as 
e — > 0, where H(a, 00) is the limit of H(a,z) when z —> 00 in any part of 
the complex plane. This limit is given by Eq. 1A91 . C±( a > bi —00) denotes 
the limit of C±(°i 6, z) when z — > 00 in the left complex half-plane. In Eq. 
1A18I . this limit is expressed in terms of (±(a, 6, +00). The integral along 
r_R tends to when R — > 00. From the theorem of residues, we have for 
any z£C\iR 

?7o,±(a,6,z) = _ — (It l)(C±)(a,6,+oo) 
VI — « 

- z {Y{-$t(z)]H(a, -z)C±(a, 6, z) =F Y{$t(z)]H(a, z)C ± (a, 6, -2)} . (27) 



G 




Figure 1: Contour for the calculation of ?7o,±( a j b, z). 



The end of the calculation of £Q follows from Eqs. lfH?t . ipty and ll2"?|l . 

Changing z into 1/z, we get 



£Q(a,6,z) = 



-0-(a, 6, +oo)m+(o, 6, 1/z), 



(28) 



a result valid in principle on C \ iR. In fact, it can easily be seen that 
the limits of this function on both sides of the imaginary axis coincide. 
This makes it possible to extend the relation 12811 by continuity on iR, 
including 0. Such an extension will not be used subsequently. 



3 Inverting the Laplace transform 

The inversion formula 1111 reads here, for any r G]0, 6[ 



1 



c+ioo 



Q(a,b,r) = — 2Q(a,b, z) exp(rz)dz (29) 

2i7T Jc — ioo 

(c / 0). Substituting 12811 in this expression and replacing u+(a,b, 1/z) 
by its definition 1201 . we obtain 

Q(a,b,r) = £-(«. b,+oo) 

Vl — a 

1 rc'+too j 
X—/ £r(o,l/ Z )C+(o,6,-l/a;)exp(Tz) — , (30) 

Z17V ./ c /_^ 00 2 



where c' = |c| > 0. The integral on the right-hand side can be calculated 
with the help of the residue theorem applied to a contour in the left half- 
plane, owing to the presence of the exponential. A difficulty arises because 
c > and the imaginary axis should not be crossed since the functions £± 
are discontinuous on this axis. Nevertheless, one can derive the identity 



H(a, l/z)Cn(a, 6, —1/z) exp(rz)- 



z 



-j H(a,l/z)(C + )-(a,b,-l/z)eMrz) T (31) 
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from the contour of Fig. 2(a), since the integral on the half-circle T £ tends 
to when e — > 0, due to Eq. 1A18I . The last integral is a Cauchy principal 
value at 0. Then the limits of the functions £± on the left-hand side of 
the imaginary axis are expressed in terms of their limits on the right-hand 
side with the help of 12511 . which shows that 



1 f +icx > dz 
— 4 H(a, l/*)(C+)-(a, b, -1/z) exp(rz)- = 

J _ iQO Z 

5S L H{ " 



l/«)(C+) + (a, b, -l/ Z )[exp(T0) + exp((6 - r)z)] 



dz 



(32) 



The last integral is calculated by displacing the line of integration on the 



(a) 



O 









(b) 






-k(a) 







> 


> 






> 











Figure 2: Contours (a) and (b) for the proof of Eqs. (pTT|) and ffify. 

left side of the imaginary axis, between the origin and the point of abscissa 
— k(a) where the function z — > H(a, 1/z) diverges: see the comments after 
Eq. 1A8I in Appendix A, which contains the definition of k(a). Applying 
the residue theorem to the contour of Fig. 2(b) and observing that the 
integral on r e tends to —H(a, oo)£+(a, b, +oo) when e — ► 0, one obtains 



i r+iao 

— j- H{a, l/z)((+) + (a, b, -l/*)[exp(rz) + exp((6 - r)z)]- 

<+(a,6,+oo) (33) 



, dz 
z 



1 



VT-a 

1 rc"+icc dz 
+ — / H(a,l/zX+(a,b,-l/z)[exp{Tz)+exp((b-T)z)]— , 

J c" — ioo % 

where c" g ] — k(a), 0[. For r £ ]0, b[, the last integral coincides with the 
function —F+(a,b,r) introduced in the Appendix B, Eq. 1BH . 
According to Eqs. 13011 - 133(1 . we finally have 

Q(a,b,r) = -== C-("A+oo)[ 2. (+(a,b, +oo) -F+(o,6,t)], 

(34) 

or, using Eq. !B3t 

Q{a,b,r) = [1 + F_(o,&,Q)][l + F+(a,6,0) - F+(a,b,r)]. (35) 
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This expression is valid a priori over ]0, 6[, but it can be extended by 
continuity over [0,6], since the functions F± are continuous on the right 
at and on the left at 6. One has in particular 



Q(a, b, 0) = Q(a, 6, 6) = 1 + F_ (a, 6, 0). 



(36) 



These results are similar to those of [5J, where details concerning the 
numerical evaluation of the function Q can be found. Ten-figure tables of 
the function (1 — a)Q are given there. The functions F± are computable 
once their definition <BH has been transformed by the method of residues, 
as outlined in Appendix B. 



Using the finite Laplace transform, we derived in a concise way the ex- 
pression 1351 of the solution to Eq. Jfjjl, which is appropriate for numerical 
evaluation 5 . Our main objective was to come directly to this expression, 
with emphasis on the manner in which the four auxiliary functions of the 
problem were generated. They are (i) the function T = T{a,z), which 
characterizes the solution to Eq. © in an infinite medium (the range 
[0,6] is replaced by R), (ii) the function H — H(a,z) which expresses its 
solution in a half-space (6 = oo), and (iii) the functions (± = (±(a,b,z) 
that complete the previous ones in the finite case (6 < oo). The main 
properties of these functions are summarized in Appendix A. We note 
that the T-function can be expressed in terms of elementary transcenden- 
tal functions [Eq. IA3I 1 . the H- function is defined explicitly by an integral 
on [0, 1] [Eq. IA10I 1. and the functions f± are defined implicitly as the 
solutions to Fredholm integral equations of the second kind [Eq. 1A13H . 
It seems that the problem JSJ has no exact closed-form solution for 6 < oo. 
The approach of this article is appropriate for solving integral equations 
of the form with a kernel defined by and any free term. It also 
applies to convolution kernels defined by functions more general than @, 
for example of the form 



where I is an interval of R and \E r a function ensuring the existence of the 
integral. This class of kernels models scattering processes more complex 
than the one considered in this article, which corresponds to I = ]0, 1] and 



*(*) = (l/2)(l/x). 

Appendix A. The auxiliary functions as- 
sociated to the problem (JBI) 

These functions are T = T(a, z), H = H(a, z), and £± = C±(°> b, z), where 

zee. 



4 Conclusion 




(37) 
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A.l The T-function 



This function is defined by 

r(a lZ ):=l-^/ +X -^- [z*±l). (Al) 

It describes the multiple scattering of photons in an unbounded medium 
for the adopted scattering law. It is defined in the whole complex plane C, 
except at ±1, provided that the integral is calculated in the sense of the 
Cauchy principal value when z G ] — 1, It is thus sectionally analytic 
in C \ [— 1, +1], its limits on both sides of the segment ] — 1, +1[ being 
given from Plemelj formulae |5j by 

T^a.u) =T(a,u) ±wr(a/2)u (-K u < +1). (A2) 

Here, T + (a,u) [resp. T~(a,u)] denotes the limit of T(a,z) when z tends 
to it G ] — 1, +1[ from above (resp. below) the real axis, and 

T(a I u) = l-^«ln(ii^) (-Ku<+1). (A3) 
2 1 — u 

The roots of the characteristic equation T(a, z) — are important for 
solving integral equations of the form Q with kernel defined by © . When 
< a < 1, this equation has four non-zero roots in C, namely two pairs of 
opposite real numbers since the T-function is even |1U|. There is a unique 
root strictly greater than unity, denoted by l/k(a) (0 < k(a) < 1), which 
is calculated by solving the transcendental equation 



T(M/*(a)) = l-y_]» 



1 + k(a) 



l-fc(a)J 



Its unique solution k(a) in ]0,1[ is given by 



k(a) = vl — a exp 

where 



1 at \ du 
o(a, u) — 



= 0. (A4) 



(A5) 



6(a,u) := -&rg[T + (a,u)] CO < u < 1). (A6) 

TV 

Here, arg(z) is the principal value of the argument of z G C. It can be 
computed using the atan2-function, since arg(i + iy) = atan2(y, x). 

A. 2 The H -function 

This function is the unique solution, analytic in the right complex half- 
plane, of the integral equation 

T(a,z)H(a,z) = l+~z f H(a,u)-^- (z^±l). (A7) 
2 Jo u — z 

This is a Cauchy integral equation on [0, 1], which defines the extension 
of H outside this interval. The function H satisfies the Wiener-Hopf 
factorization relation 

T(a, z)H(a, z)H(a, -z) = 1 (z ^ ±1), (A8) 
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first proved regardless of <A7t (see, e.g., |10|). then as a result of IA7t |12|. 
We note that it implies the divergence of the H-function at the solutions of 
the equation T(a, z) — located in the left complex half-plane, especially 
at z = —l/k(a). The relations 1A7I - 1A8I also mean that the H- function 
is bounded at infinity, with limit 

ff(a,oo) = ^=L= (A9) 
V 1 — a 

whatever the region of the complex plane where z goes to infinity. This 
results from the fact that T(a, oo) = 1 — a. 

Many analytical expressions of the function H have been derived since 
the introduction of this function in the thirties. The following expression 
is valid in the right complex half-plane [5j 



k{a) l + z 

H(a,z) = exp 

Vl - a 1 + k(a) z 



-z / 9(a,u) 
Jo 



du 



u + z 



(»(*) > 0), 



(A10) 

where 6 is the function defined by 1A6I . We used it to compute H(a,z) 
over [0, 1] and at l/k(a), a necessary step for computing the functions F± 
appearing in our final result 13511 : see the Appendix B. 

A. 3 The functions C± 

These functions are sectionally analytic in C \ iR, with limits (C±) + an d 
(C±)~ on the right- and left-hand sides of the imaginary axis, respectively. 
They are denned by the relations 



C±(a,b,z) = 1 ± — 



H(a,~l/z') 
2i7v J_ loo H(a, l/z>) 



xe X p(-bz')(( ± ) + (a,b,l/z')— -,(A11) 

1+zz 

which can be transformed into Fredholm integral equations of the second 
kind on [0, 1] using the residue theorem |13|. We suppose here the exis- 
tence and uniqueness of the solution to Eq. lAlH . and refer to |13| for 
the proof of the relation 

i[C+(a, b, z)C-(a, b, -z) + C-(a, b, z)C+(a, b, -*)] = 1, (A12) 

which is valid over C \ iR. 

The transformation of Eq. lAlU into a couple of Fredholm integral equa- 
tions on [0, 1] is described in the section V of We reproduce the 
relation (71) of this article, valid for any z in C \ {] — 1, 0[U { — l/k(a)}} 

2kz 

C±(a,b,z) = l±M±(o,&) : " 



1 + kz 



ii 



We have introduced in the right-hand side 

M ± (a,b) :=«(o,6)C±(o,6,l/fc(o)), (A14) 

where 

qia ^ : =l H ^S(a)) eM - k{a)b) ' (A15) 

R(n) ._ k ( a ) _ l-fc 2 (a) , A1fl x 
K(a) - T'(a,l/fc(a)) ~ fe2(a)+a-l' 1 &j 

and the function g in the integral term is 

: = TfTt Lrj 7 = 7pr ( \~T7 7w (0< u < !)• ( A1 ?) 

l + (a,u)l (a,u) 1 z (a,u) + (irau/ iy 

Let z go to infinity in both members of Eq. HA13I , first in the left complex 
half-plane, then in the right one. Writing C±(a, b, — oo) and C±( a > b, +oo) 
the limits obtained in this way, we have 

C±(a, b, -oo) = (1 =f l)C±(a, b, +oo). (A18) 

Consequently 

C-(a,6, +oo)C+(a, b, +oo) = 1, (A19) 

due to JXT2J. 

The numerical calculation of the functions £± is carried out starting from 
Eq. 1A13L which consists of two Fredholm integral equations on [0,1]. 
A difficulty arises because the coefficients M± (a, b) contain the unknown 
values of £± at l/k(a). It is overcome by putting z = l/k(a) into Eq. 
IA13t to calculate M±(a, b), which is next substituted into the right-hand 
side of llATSTl . The expression of the functions £± in the right half-plane 
becomes |13| 

(±(a,b,z) = p±(a,b,z) ±M±(a,b)a±(a,b,z), (*R(z) > 0), (A20) 

where p± and a± are solution to the integral equations 

a f 1 q(a,u) , , , . . . du , . „ H . 

~ ' yv ' ; exp(-b/u)p±(a,b,u) , (A21) 



ll) u + z 



/ , s 2k(a)z , a o(a, u) , , , , , , . du 

a ± (a,b,z) = )/ ±-z / " v ' exp{-b u)a±(a,b,u) 
1 + k(a)z 2 J n^(a,u) 



u + z 
(A22) 

Once they have been solved on [0, 1], these equations define the functions 
z — > p±(a, b, z) and z — > <r±(a, fe, z) everywhere in the right complex half- 
plane, in particular at z = l/k(a). Whence the coefficients M±(a,6) from 
Eqs. <ATH and llA"20ll arc 

±V ' ; l + q(a,b)a ± {a,b,l/k{a))' y ' 

q(a,b) being defined by (|A15^ - (|A16fl . 

Calculating the functions £± thus requires the numerical solution of Eqs. 
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1A21I - IA22I . In fact, it can be proved that the solution of Eqs. 1A22II 
comes down to that of Eqs. 1A211 . This means that the functions a± are 
analytically expressible in terms of the functions p±. It follows that only 
Eqs. 1A211 have to be solved numerically in view of computing the func- 
tions £±. The details are not given here, since there is no inconvenience, 
at the level of accuracy, to solve both sets of equations numerically. This 
results from the regularity of the free terms and kernels of Eqs. 1A21I - 
1A22I on [0, 1] and [0, 1] x [0, 1] respectively. Their solution p± and a± 
are easily computable since they are regular and smooth on [0, 1], i.e., not 
only continuous, but differentiable everywhere in [0, 1], including at 0. We 
plan to publish ten-figure tables of these functions in the future. 

Appendix B. The functions F± = F±(a, b, r) 

These functions are defined by the relation 

F ± (a,b,r) ~ ~-[5(r)±5(b-r)]- — / H(a, l/«)C±(o, b, -1/z) 
2 2itt J c _ ioo 

x [exp(rz) ± exp((b - t)z)}— (Bl) 

where —k(a) < c < 0. <5(t) and S(b — r) are the Dirac distributions at 
and b, respectively. 

B.l Three properties of the functions F± 

We first prove the following relations: 

F±(a,b,b-r) = ±F ± (a,b,r), (B2) 



F ± (a,b,Q) = -^=C±{a,b,+oo)-l, (B3) 

(1- a) [l + F+(a,b,Q)] [1 + F_ (a, &, 0)] = 1. (B4) 

Equation 1B2I is evident and Eq. 1B4II is a direct consequence of Eqs. 
(iBlt and JXTSJi. To derive we note that 

F±(o,6,0)=---— / H(a,l/z)C±(a,b,-l/z)[l±e W (bz)]— . 



(B5) 

The integral is transformed with the help of the residue theorem applied 
to the contour of Fig. 3, within which the integrand is analytic. Since the 
integral along the half-circle F e tends to (l/2)H(a, oo)£±(a, b + oo)(l ± 1) 
when e — > 0, we obtain 

F ± (a, b, 0) = -1 + i(l ± l)H(a, oo)C±(o, 6, +co) 

- ^- H( - a > V^)(C±) + («, b, -l/z)[l ± exp(6z)] — . (B6) 
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Figure 3: Contour for the calculation of F±(a, b, 0). 

Put 

I f+ioa j 

(«)=-— y i?(a,l/ 2 )(C±) + (a,6,-l/z) T , 

I r+ica j 

(/3)=T^ y #(a,l/ 2 )(C±) + (a,b,-l/z)exp(fe)- 

so that the integral term in HB6t is (a) + (/3). From Eq. 125H with zq = 1/z, 
we have 

1 f+i°° 1 

(0) = -T~ H ( a ' - V^)(C±) + (a, b, 1/*)- 

1 /' +i00 rfz 
+ ^-/ ff(a,-l/r)(C±)-(a,6 1 l/«)-, 

2 ™ J-ioo Z 

= +2^y H(a,l/^)(C ± ) + (a,& ! -l/,) T 

+ ff(a,-l/ Z )(C±)-(a ) 6,l/«)^. 

Therefore 

(«) + (/3) = +2^ / ff(a,-l/ Z )(C±)-(a,6,l/ 2! )^. 

Now integrate the function z — > -W(a, — l/z)£±(a, fa, l/z)/z along the con- 
tour of Fig. 4. When e — > 0, the integral over the half-circle r e tends to 
-(l/2)J J(q, oo )C±(a i b, -co) = -(l/2)i7(a,oo)(l=F l)C±(a,6,+oo) owing 
to Eq. HA18L and the integral over F_r tends to 1/2 when 7? goes to 
infinity. Hence 

(a) + 09) = i(lTl)-ff(o,c»)C±(a,6, +oo) - i 
and Eq. |B6) leads to {B^ since H(a,oo) = 1/yi - a from Eq. l|X9t . 
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Figure 4: Contour for the calculation of (a) + (/3). 
B.2 Numerical calculation of F±(a, b, r) 

It is carried out by replacing H(a, 1/z) by l/[T(a, l/z)H(a, — 1/z)] in the 
integral of IB1I . according to 1A8L This reveals the pole z — —k(a) of the 
integrand in the left complex half-plane, together with the cut ] — oo, —1] 
insuring that the function z — > T(a, 1/z) is single- valued in this half- plane. 
This function is then sectionally analytic in the complex plane cut along 
] — oo, — 1[U] + 1, +oo[, with limits T ± (a, 1/u) on the cut given by 1A2H . 
Next, the residue theorem is applied to the contour of Fig. 5, which passes 
round the cut. The unique pole within the contour is z — —k[a), with 



[R(a)/H(a, l/k(a))]C±(a,b, l/fc(a))[exp(-fc(a)r) ± exp(-fc(a)(6 - r))], 




Figure 5: Contour for the calculation of the functions F±. 



residue 



15 



since R(a) is defined by 1A16II . The integral on the half-circle Fr tends 
to [S(t) ± S(b — r)]/2 when R — > oo, the integral on T_i tends to with 
the radius, and the limit of the integral over AB+CD is 

2 J 1 H(a l/J) C±(a ' ' V«)[exp(-T«) ± exp(-(6 - r)t>)]^, 

where g is defined by Eq. (IA17II . Putting it = 1/v in the above integral, 
one obtains 

F± fa^r)= g( ^ (a)) C fe (a,M/*(a)) 

X [exp(— fc(a)r) ± exp(— k(a)(b — r))] 

a /• 1 |(^ c±(aj6jU)[exp( _ r/u)±exp( _ (fe _ r)/M)] ^ (B7) 
2 Jo H(a,u) 

For given functions fe(a), (a, u) and C,± (a, b, u), the numerical evaluation 
of the functions F± by means of this expression is easy, except when o — * 1 
as regards the integrated term, and when r — > + or b~ as regards the 
integral term. These difficulties are overcome in [Sj. 
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